In [ ]:
import os
from os import path
# Third-party
from astropy.io import ascii, fits
import astropy.coordinates as coord
from astropy.time import Time
import astropy.units as u
import matplotlib as mpl
from matplotlib.gridspec import GridSpec
import matplotlib.pyplot as plt
import numpy as np
%matplotlib inline
from scipy.stats import beta
import h5py
from thejoker import JokerSamples
from thejoker.sampler import JokerParams, TheJoker
from thejoker.plot import plot_rv_curves
from twoface.config import TWOFACE_CACHE_PATH
from twoface.db import (db_connect, AllStar, AllVisit, AllVisitToAllStar,
StarResult, Status, JokerRun, initialize_db)
from twoface.data import APOGEERVData
from twoface.plot import plot_data_orbits, plot_two_panel
from twoface.mass import m2_func
from twobody import KeplerOrbit
In [ ]:
Session, _ = db_connect(path.join(TWOFACE_CACHE_PATH, 'apogee.sqlite'))
session = Session()
samples_file = path.join(TWOFACE_CACHE_PATH, 'apogee-jitter.hdf5')
mcmc_samples_file = path.join(TWOFACE_CACHE_PATH, 'apogee-jitter-mcmc.hdf5')
In [ ]:
troup = ascii.read('../../../papers/thejoker-paper/data/troup16-dr12.csv', format='commented_header')
In [ ]:
figures_path = '../../paper/1-catalog/figures/'
In [ ]:
in_ = []
not_in = []
with h5py.File(samples_file) as f:
for apogee_id in troup['APOGEE_ID']:
if apogee_id in f:
in_.append(apogee_id)
else:
not_in.append(apogee_id)
len(in_), len(not_in)
Q: Why aren't all Troup stars in our sample?
A: logg cut
In [ ]:
# data = star.apogeervdata()
# _ = data.plot()
In [ ]:
# allvisit = fits.getdata('/Users/adrian/data/APOGEE_DR14/allVisit-l31c.2.fits')
In [ ]:
# visits = allvisit[allvisit['APOGEE_ID'] == star.apogee_id]
# visits['STARFLAGS']
In [ ]:
def troup_to_orbit(row):
P = row['PERIOD']*u.day
e = row['ECC']
K = row['SEMIAMP'] * u.m/u.s
a_K = P * K / (2*np.pi) * np.sqrt(1 - e**2)
orbit = KeplerOrbit(P=P, e=e, a=a_K,
omega=row['OMEGA']*u.rad,
t0=Time(row['T0'], format='jd'),
i=90*u.deg, Omega=0*u.deg, M0=0*u.deg)
orbit._v0 = (row['V0'] + row['SLOPE']*row['T0']) * u.m/u.s
return orbit
In [ ]:
plot_path = '../../plots/troup-compare/'
if not path.exists(plot_path):
os.makedirs(plot_path)
In [ ]:
with h5py.File(samples_file) as f:
for i, apogee_id in enumerate(in_):
samples = JokerSamples.from_hdf5(f[apogee_id])
star = AllStar.get_apogee_id(session, apogee_id)
data = star.apogeervdata(clean=True)
troup_row = troup[troup['APOGEE_ID'] == apogee_id]
troup_orb = troup_to_orbit(troup_row)
fig = plot_two_panel(star, samples, title=star.apogee_id,
plot_data_orbits_kw=dict(highlight_P_extrema=False,
n_times=16384,
plot_kwargs=dict(color='#666666')))
ax1, ax2 = fig.axes
# over-plot Troup orbit
t2 = Time(np.linspace(*ax1.get_xlim(), 10000), format='mjd')
ax1.plot(t2.tcb.mjd, troup_orb.radial_velocity(t2).to(u.km/u.s),
marker='', color='tab:orange', alpha=0.5)
ax2.scatter(troup_row['PERIOD'], troup_row['ECC'],
marker='+', linewidth=2., s=100, color='tab:orange',
label='Troup')
fig.savefig(path.join(plot_path, '{0}.pdf'.format(apogee_id)))
plt.close(fig)
In [ ]:
classes = dict()
classes['unimodal'] = ['2M04411627+5855354', '2M19405532+2401157', '2M03080601+7950502', '2M19134104-0712053']
classes['multimodal'] = ['2M00295684+6356284', '2M19453527+2333077', '2M19114515-0725486', '2M19292561+2626538']
classes['data changed'] = ['2M18591837-0401083', '2M19105197+2845422']
In [ ]:
rc = {
'axes.labelsize': 18,
'xtick.labelsize': 14,
'ytick.labelsize': 14
}
with mpl.rc_context(rc):
with h5py.File(samples_file) as f, h5py.File(mcmc_samples_file) as mcmc_f:
for label, ids in classes.items():
# fig, axes = plt.subplots(len(ids), 2, figsize=(12, 4*len(ids)))
fig = plt.figure(figsize=(8, 2.375*len(ids)))
gs = GridSpec(len(ids), 3)
for i, apogee_id in enumerate(ids):
if apogee_id in mcmc_f:
samples = JokerSamples.from_hdf5(mcmc_f[apogee_id])
print('mcmc')
else:
samples = JokerSamples.from_hdf5(f[apogee_id])
print('thejoker')
star = AllStar.get_apogee_id(session, apogee_id)
data = star.apogeervdata(clean=False)
samples.t0 = data.t0
troup_row = troup[troup['APOGEE_ID'] == apogee_id]
troup_orb = troup_to_orbit(troup_row)
ax1 = fig.add_subplot(gs[i, :2])
ax2 = fig.add_subplot(gs[i, 2])
axes = [ax1, ax2]
scatter_kw = dict()
if i == 0:
scatter_kw['label'] = 'The Joker'
plot_two_panel(data, samples, axes=axes,
plot_data_orbits_kw=dict(highlight_P_extrema=False,
n_times=16384,
relative_to_t0=True,
plot_kwargs=dict(color='tab:blue')),
scatter_kw=scatter_kw)
ax1.set_xlabel('')
ax2.set_xlabel('')
# over-plot Troup orbit
_t0 = data.t.tcb.mjd.min()
t2 = Time(np.linspace(*ax1.get_xlim(), 10000) + _t0, format='mjd')
ax1.plot(t2.tcb.mjd - _t0, troup_orb.radial_velocity(t2).to(u.km/u.s),
marker='', color='tab:orange', alpha=0.5, zorder=-10, linewidth=2)
ax2.scatter(troup_row['PERIOD'], troup_row['ECC'],
marker='+', linewidth=2., s=100, color='tab:orange',
label='Troup+2016 fit', zorder=-10, alpha=0.75)
if i == 0:
ax2.legend(loc='upper left', fontsize=10)
# add text
xlim = ax1.get_xlim()
ylim = ax1.get_ylim()
ax1.text(xlim[0] + (xlim[1]-xlim[0])/20,
ylim[1] - (ylim[1]-ylim[0])/20,
apogee_id, fontsize=12, ha='left', va='top')
ax1.set_xlabel(r'${\rm BMJD} - t_0$ [day]')
ax2.set_xlabel('period, $P$ [day]')
fig.tight_layout()
fig.savefig(path.join(figures_path, 'troup-{0}.pdf'.format('-'.join(label.split()))))
plt.close(fig)
In [ ]:
fig, ax = plt.subplots(1, 1, figsize=(8, 5))
sub = troup[troup['SNR']>100]
cb_in = ax.scatter(sub['PERIOD'], sub['ECC'], marker='.', cmap='magma_r',
c=sub['NVISITS'], vmin=3, vmax=20)
ax.set_xscale('log')
ax.set_xlim(1, 2000)
ax.set_xlabel('$P$ [{0:latex_inline}]'.format(u.day))
ax.set_ylabel('$e$')
ax.set_title('SNR > 100 - Troup')
cb = fig.colorbar(cb_in)
cb.set_label('$N$ visits')
In [ ]:
sub = troup[(troup['SNR'] > 100) & (troup['PERIOD'] > 8)]
bins = np.linspace(0, 1, 13)
fig, axes = plt.subplots(1, 2, figsize=(10, 5), sharex=True)
mask = sub['PERIOD'] < 20
axes[0].hist(sub[mask]['ECC'], bins=bins, normed=True, alpha=0.8);
axes[0].set_title(r'$8 < P < 20\,{{\rm d}}$ ({0} stars)'.format(mask.sum()))
axes[1].hist(sub[~mask]['ECC'], bins=bins, normed=True, alpha=0.8);
axes[1].set_title(r'$P > 20\,{{\rm d}}$ ({0} stars)'.format(np.logical_not(mask).sum()))
ecc = np.linspace(0, 1, 100)
for ax in axes:
ax.plot(ecc, beta.pdf(ecc, 0.867, 3.03), marker='', label='prior')
ax.set_xlabel('eccentricity, $e$')
fig.suptitle('Troup', y=1.02, fontsize=20)
fig.tight_layout()
In [ ]:
fig, ax = plt.subplots(1, 1, sharex=True, figsize=(5,4))
n, bins, _ = ax.hist(troup['OMEGA'], bins='auto');
binc = (bins[:-1]+bins[1:])/2.
ax.errorbar(binc, n, np.sqrt(n), marker='', linestyle='none')
ax.set_xlabel('$\omega$ [rad]')
In [ ]:
fig, ax = plt.subplots(1, 1, figsize=(6, 5))
sub = troup[ (troup['SNR'] > 100) &
(troup['FE_H'] > -999)]
print(len(sub))
ax.errorbar(sub['PERIOD'], sub['FE_H'], yerr=sub['FE_H_ERR'],
linestyle='none', marker='.', color='k')
ax.set_xscale('log')
ax.set_xlim(1, 2000)
# ax.set_xlabel('[Fe/H]')
# ax.set_ylabel(r'$M_{2, {\rm min}}$ ' + '[{0:latex_inline}]'.format(u.Msun))
# ax.set_title('log$g$ < 3.25, $\chi^2$ < 30')
fig.tight_layout()
In [ ]:
In [ ]: